#ifndef MCMC_SPARSEMATRIX_H
#define MCMC_SPARSEMATRIX_H

#include "mcRaster.h"

template<class Tp>
class SparseMatrix : public Raster<Tp>
{
public:
	// construction & destruction
	SparseMatrix();
	~SparseMatrix();
	void FreeSparse(void);

	// operator
	SparseMatrix<Tp>& operator=(const SparseMatrix &a);

	void	SetNumSparse(const int num);
	inline	int NumSparse(void) const {return num_sparse;};
	bool	SetValue(const int index, const Tp value, const int row, const int col);
	bool	GetValue(const int index, Tp &value, int &row, int &col) const;
	void	FillMatrix(void);
	mcmcMatrix<Tp> SparseConv(mcmcMatrix<Tp> &input_matrix);
	void	Transfer(const mcmcMatrix<Tp> &a);
	void	Flip(void);

protected:
	Tp		*psparse_data;		// valid real data
	int		*prow;				// valid cell in row
	int		*pcol;				// valid cell in column
	int		num_sparse;			// number of sparse cell
};


// construction & destruction
template<class Tp>
SparseMatrix<Tp>::SparseMatrix()
{
	psparse_data	= NULL;
	prow			= NULL;
	pcol			= NULL;
	num_sparse		= 0;
}

template<class Tp>
SparseMatrix<Tp>::~SparseMatrix()
{
	FreeSparse();
}

template<class Tp>
void SparseMatrix<Tp>::FreeSparse(void)
{
	if (psparse_data != NULL)	delete []psparse_data;
	psparse_data = NULL;

	if (prow != NULL)	delete []prow;
	prow = NULL;

	if (pcol != NULL) delete []pcol;
	pcol = NULL;

	num_sparse = 0;
}


// operator
template<class Tp>
SparseMatrix<Tp>& SparseMatrix<Tp>::operator=(const SparseMatrix &a)
{
	if (this != &a)
	{
		int row,col;
		Tp	value;

		SetNumSparse(a.NumSparse());
		for (int i=0; i<NumSparse(); i++)
		{
			value = a.GetValue(i,row,col);
			SetValue(value,i,row,col);
		}

		// copy the matrix itself
		SetDimension(a.rows(),a.cols());
		for (int j=0; j<rows(); j++)
			for (i=0; i<cols(); i++)
				(*((mcmcMatrix<Tp>*) this))(j,i) = a.Data(j,i);
	}

	return *this;
}

// operations
template<class Tp>
void SparseMatrix<Tp>::SetNumSparse(const int num)
{
	if (num>0 && num!=NumSparse())
	{
		FreeSparse();

		psparse_data = new Tp[num];
		prow		 = new int[num];
		pcol		 = new int[num];

		num_sparse = num;
	}
}

template<class Tp>
bool SparseMatrix<Tp>::SetValue(const int index, const Tp value, const int row, const int col)
{
	if (index>=0 && index<NumSparse())
	{
		psparse_data[index] = value;
		prow[index]			= row;
		pcol[index]			= col;

		return true;
	}
	return false;
}

template<class Tp>
bool SparseMatrix<Tp>::GetValue(const int index, Tp &value, int &row, int &col) const
{
	if (index>=0 && index<NumSparse())
	{
		value = psparse_data[index];
		row = prow[index];
		col = pcol[index];

		return true;
	}
	return false;
}

template<class Tp>
void SparseMatrix<Tp>::FillMatrix(void)
{
	Tp	value;
	int row,col;

	InitValue(0);
	for (int i=0; i<NumSparse(); i++)
	{
		GetValue(i,value,row,col);
		(*((mcmcMatrix<Tp>*) this))(row,col) = value;
	}
}

template<class Tp>
mcmcMatrix<Tp> SparseMatrix<Tp>::SparseConv(mcmcMatrix<Tp> &input_matrix)
{
	mcmcMatrix<Tp>	I;
	int		i,j,i1,j1,i2,j2,m1,n1,m2,n2,is;
	Tp			value;

	I.SetDimension(input_matrix.rows(), input_matrix.cols());
	I.Zero();

	m1 = rows();
	n1 = cols();

	m2 = input_matrix.rows();
	n2 = input_matrix.cols();

	for (j=0; j<m2; j++)
		for (i=0; i<n2; i++)
		{
			for (is=0; is<num_sparse; is++)
			{
				j1 = prow[is];
				i1 = pcol[is];
				value = psparse_data[is];

				j2 = j-(j1-(m1-1)/2);
				i2 = i-(i1-(n1-1)/2);
				if (j2>=0 && j2<m2 &&i2>=0 && i2<n2)
					I(j,i) = I(j,i)+input_matrix.Data(j2,i2)*value;
			}
		}
	return I;
}

// transfer a normal matrix to this sparse matrix
template<class Tp>
void SparseMatrix<Tp>::Transfer(const mcmcMatrix<Tp> &a)
{
	// count the number of non-zero elements
	int i,j,num=0,k;

	for (j=0; j<a.rows(); j++)
		for (i=0; i<a.cols(); i++)
			if (a.Data(j,i) != 0)
				num++;

	// set this sparse matrix
	SetDimension(a.rows(),a.cols());
	SetNumSparse(num);
	k = 0;
	for (j=0; j<a.rows(); j++)
		for (i=0; i<a.cols(); i++)
		{
			if (a.Data(j,i) != 0)
				SetValue(k++,a.Data(j,i),j,i);
		}

	FillMatrix();
}

// flip this matrix
template<class Tp>
void SparseMatrix<Tp>::Flip(void)
{
	Tp	value;
	int row,col;

	for (int i=0; i<NumSparse(); i++)
	{
		GetValue(i,value,row,col);
		SetValue(i,value,rows()-1-row,cols()-1-col);
	}

	FillMatrix();
}

#endif
